#ifndef _HELPER_H_
#define _HELPER_H_

#include "Spline.h"
#include "SplineHelper.h"
#include <iostream>
#include <vector>
#include <utility>
#include <cmath>


typedef double (*fpd)(double);

double Bfunc(const double x){
     return 1.0/(1+25*x*x);
}

auto DataGenerator(fpd f, const double start, const double end, const
int num){
     std::vector<double> v1,v2;
     const double h = (end - start)/num;
     double t = start;
     std::vector<Vec<double,2> > res;
     for (int i = 0 ; i <= num ; i++){
       res.push_back({t,f(t)});
       t += h;
     }
     return res;
}

double Berror(const Spline<1,4,ppForm>& sp, const int N){
  const double h = 2.0/N;
  double t = -1.0 + h/2;
  double Max = 0;
  for (int i = 0 ; i < N ; i++){
    Max = std::max(Max,fabs(sp(t)[0]-Bfunc(t)));
    t += h;
  }
  return Max;
}

double Cfunc(const double x){
  return 1.0/(1+x*x);
}

auto Cgenerator1(){
  std::vector<double> x;
  std::vector<double> y;
  x.push_back(-5);
  y.push_back(5.0/338);
  for (int i = -5 ; i <= 5 ; i++){
    x.push_back(i);
    y.push_back(Cfunc(i));
  }
  x.push_back(5);
  y.push_back(-5.0/338);
  return std::make_pair(x,y);  
}

auto Cgenerator2(){
  std::vector<double> x;
  std::vector<double> y;
  x.push_back(-5);
  y.push_back(Cfunc(-5));
  for (int i = -5 ; i < 5 ; i++){
    x.push_back(i+0.5);
    y.push_back(Cfunc(i+0.5));
  }
  x.push_back(5);
  y.push_back(Cfunc(5));
  return std::make_pair(x,y);  
}


template <int Ord>
double Derror(const Spline<1,Ord,BForm>& Bs, const double t){
  return fabs(Bs(t) - Cfunc(t));
}

double Efunc1(const double x){
  double tmp = std::sqrt(fabs(x))-std::sqrt(3-x*x);
  return 2.0*tmp/3;
}

double Efunc2(const double x){
  double tmp = std::sqrt(fabs(x))+std::sqrt(3-x*x);
  return 2.0*tmp/3;
}

auto Egenerator(const int num){
  assert(num%2 == 0);
  std::vector<Vec<double,2> > vdata1;
  std::vector<Vec<double,2> > vdata2;
  vdata1.push_back({0,-2.0*std::sqrt(3)/3});
  vdata2.push_back({0,-2.0*std::sqrt(3)/3});
  const int N2 = (num - 4)/3;
  const int N1 = num/2 - 2 -N2;
  const double h1 = std::sqrt(3)/(N1+1);
  const double h2 = std::sqrt(3)/(N2+1);
  double t = 0;
  for (int i = 0; i < N1 ; i++){
    t += h1;
    double tmp = Efunc1(t);
    vdata1.push_back({t,tmp});
    vdata2.push_back({-t,tmp});
  }
  vdata1.push_back({std::sqrt(3),2.0*std::sqrt(std::sqrt(3))/3});
  vdata2.push_back({-std::sqrt(3),2.0*std::sqrt(std::sqrt(3))/3});
  t = std::sqrt(3);
  for (int i = 0; i < N2 ; i++){
    t -= h2;
    double tmp = Efunc2(t);
    vdata1.push_back({t,tmp});
    vdata2.push_back({-t,tmp});
  }
  vdata1.push_back({0,2.0*std::sqrt(3)/3});
  vdata2.push_back({0,2.0*std::sqrt(3)/3});
  return std::make_pair(vdata1,vdata2);
}

#endif // _HELPER_H_
